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Abstract 


This book presents a collection of efficient machine language 
programs for advanced applications with the TMS320C30 family 
of digital processors. These programs include the following 
categories: 


Normal precision floating point math functions 
Extended precision floating point math functions 
Integer arithmetic routines 

Vector utility routines 


Radix 2 FFT routines 


Oooo vO DO 


Linear algebra routines 


The names and short descriptions of these routine categories (and 
special notations) are included. 


The book contains detailed information about: 
Extended vs. Normal Precision 

Program utilization 

Function approximation techniques 

Math function details 

Integer arithmetic program details 


Vector utility routines 


Ooo O O OD DO 


Ftp routines 
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Q Linear algebra routines 


The book concludes with a list of references and an appendix of 
source code for the described routines. 
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Introduction 


This report presents a collection of efficient machine language programs for advanced 
applications with the TMS320C30. These programs provide basic math and transcenden- 
tal functions. Other routines include vector functions, FFTs and linear algebra. 


Library Overview 
The set of programs fall into six categories: 


I. Normal precision floating point math functions, 
II. Extended precision floating point math functions, 
Il. Integer arithmetic routines, 

IV. Vector utility routines, 

V. Radix 2 FFT routines, and 

VI. Linear algebra routines. 


Categories I and II are programs which implement a minimal set of elementary 
mathematical functions for advanced applications. In these categories, the functions FPINV 
and SQRT are improved versions of the programs in the TMS320C3x User’s Guide (1). 
In category III, IMULT and IDIV are improved versions of the programs EXTMPY and 
DIVI in [1]. In category IV, *FMIEEE and *TOIEE are array versions of the TOIKEE 
and FMIEEE scalar programs from the User’s Guide. 


The names and short descriptions of these routines use some special notation: 


Categories I and II: xd — indicates that the relative accuracy of the im- 
plemented function is x decimal digits. 
Categories IV and VI: * — program name prefix stands for M or R. 


M — selects the memory based parameter entry point. 

R — selects the register based parameter entry point. 

Categories II and VI: X — indicates the extended precision program 
version. 


Consult the program source listings for more details. 
The following are brief descriptions of the programs by category: 
I. Normal floating-point (32-bit) math functions ($;MATH.ASM): 


A. SIN —computes a 7d sine(x) for all x in radians. 
B. Cos —computes a 7d cosine(x) for all x in radians. 
C. EXP —computes a 7d exp(x) for all |x| < 88. 
D. LN —computes a 7d In(x) for all x > 0. 
E. ATAN —computes a 7d atan(x) in radians for all x. 
F. SQRT —computes an 8d sqrt(x) for all x > 0. 
G. FPINV —computes an 8d 1/x for all x # 0. 
H. FDIV —computes an 8d x/y for all x and all y # 0. 

II. Extended-precision, floating-point (40-bit) math functions (;]MATHX.ASM): 
A. SINX —computes a 9d sine(x) for all x in radians. 
B. COSX —computes a 9d cosine(x) for all x in radians. 
C. EXPX —computes a 9d exp(x) for all |x| < 88. 
D. LNX —computes an 8d In(x) for all x > 0. 
E. ATANX —computes an 8d atan(x) in radians for all x. 
F. SQRTX —computes a 10d sqrt(x) for all x > 0. 
G. FPINVX —computes a 10d 1/x for all x # 0. 
H. FDIVX —computes a 10d x/y for all x and all y # 0. 
I. FMULTX —computes a 10d x+y for all x and y. 

III. Integer (32-bit) math routines ($;MATHI.ASM): 
A. ILOG2 —computes m = log2(n), n < 2™ for use with radix 

2 FFT programs. 
B. IMULT —computes 64-bit product of two 32-bit numbers. 
C. IDIV —computes quotient and remainder of two 32-bit 
numbers. 


IV. Vector utilities (SVECTOR.ASM): 


7 *CORMULT -— in-place computation of the complex vector pro- 
duct of two complex arrays using the complex con- 
jugate of the second array. 

B. *CONMULT -— in-place computation of the complex vector pro- 
duct of two complex arrays. 


C. *CBITREV —in-place bit reverse permutation on a complex ar- 
ray with separate real and imaginary arrays. 
D. *FMIEEE —in-place fast conversion of an IEEE array to a 


TMS320C30 array. 


E. *TOIEEE —in-place fast conversion of a TMS320C30 array to 
an IEEE array. 


F. *VECMULT — in-place multiplies a constant times an array. 
G. *CONMOV = —moves (fills) a constant into an array. 
H. *VECMOV —moves (copies) an array into another array. 
V. Radix 2 FFT routines ($FFT2.ASM): 
A. CFFFT2 —Complex DIF forward radix 2 FFT using separate 
real and imaginary arrays and 3/4 cycle sine table. 
B. CIFFT2 —Complex DIT inverse radix 2 FFT using separate 


real and imaginary arrays and 3/4 cycle sine table 
(does not include the 1/N scale factor). 


VI. Linear algebra routines ($LINALG.ASM): 


A. *SOLUIN —Solves a well conditioned system of linear equa- 
tions with any number of dependent variable sets. 
Uses no (diagonal) pivoting with normal-precision 
floating-point math. 

B. *SOLUTNX  -—Solves a well conditioned system of linear equa- 
tions with any number of dependent variable sets. 
Uses no (diagonal) pivoting with extended- 
precision floating-point math. 


Extended vs. Normal Precision 


Categories I, II, and VI represent a dual collection of programs implemented with 
32-bit single- or normal-precision TMS320C30 floating-point arithmetic, and with 40-bit 
extended-precision TMS320C30 floating-point arithmetic. Some of the normal-precision 
programs (category I, for example) have been written using the TMS320C30 RND in- 
struction for rounding to obtain the optimal precision from the standard floating point 
TMS320C30 instruction set. This has been done with a slight loss of speed. Such round- 
ing can be carefully eliminated by the user if the additional speed is necessary at the ex- 
pense of some accuracy. 


Extended-precision was implemented on the TMS320C30 by the simple implemen- 
tation of the 40-by-40 floating-point multiply routine, FMULTX. This was necessary since 
the TMS320C30 has 40-bit addition and subtraction instructions, but the multiply operates 
only on 32-bit inputs. By using the native add and subtract FMULTX and the extended- 
precision registers RO to R7, 40-bit floating-point math was effected. All 40-bit constants 
are stored in two consecutive words in memory. The first word is the normal truncated 
32-bit floating-point number. The least significant byte of the second word contains the 
remaining bottom 8 bits of the extended mantissa. The programs are coded to properly 
load extended-precision registers with these double-word constants. 


The extended-precision versions of the programs in this report may be slower than 
their normal precision counterparts. When using extended-precision results in RO from 
category II programs, note that the results may be stored in memory with or without round- 
ing, A more accurate normal-precision result will generally be obtained by rounding. You 
should never round before using an extended-precision result as input to another extended- 
precision program unless special circumstances exist. Note that truncation, not rounding, 
will occur if an extended-precision register is moved to any 32-bit register or any memory 
location. This will generally cause loss of accuracy in the amount of the value of the least 
significant bit of the mantissa. 


Program Utilization 


Since all programs in this collection are intended to be invoked by a CALL instruc- 
tion, you must have the stack pointer (SP register) appropriately set to an available memory 
area, preferably in internal RAM. Programs in categories I and II save and restore the 
data page register DP by using the stack area pointed to by SP. Programs in category 
III do not alter or use the DP register at all. The programs in categories IV through VI 
alter but do not restore the DP register. 


All of the programs in categories I through III, except for ILOG2, are implemented 
as straight line code. You may wish to disable the instruction cache while these programs 
are being executing. This will cause no loss of execution speed and will avoid flushing 
out potentially reusable instructions in the cache. It is beneficial to have the cache enabled 
when using most of the remaining programs (categories IV through VI) as they generally 
contain multi-instruction loops. 


Programs in categories IV through VI allow input through externally defined variables 
addresses. The .global references indicate these addresses, where the input variable values 
and/or addresses are located. The starting address of these memory locations is given by 
the external variable $PARAMS. All of the addresses are assumed to be in the same 
TMS320C30 memory page as $PARAMS. If this is not the case, the addresses or the 
programs should be changed assure that the DP register gets set properly. 


Programs in categories IV and VI also allow the use of registers to hold input 
parameters. The exact registers to be used are found in the program source listings. When 
using the register input entry point, refer to the program using the R prefix on the pro- 
gram name, e.g. RSOLUTN. The memory based parameter input entry uses the M prefix, 
e.g. MSOLUTN. The .global references to the R prefix entry points may be deleted if 
they are not needed. 


Function Approximation Techniques 


Categories I and II are made up of a collection of elementary mathematical func- 
tions numerically approximated using two basic methods. The functions SIN, COS, EXP, 
LN, and ATAN are approximated by using polynomials fitted to the various functions 
over a limited range of the independent variable. The functions SQRT and FPINV are 
approximated by iteratively solving a particular non-linear equation. The extended preci- 
sion versions of these programs (category II) use the same approach with extended-precision 
arithmetic and resort to more accurate polynomials or more iterations to achieve the desired 
precision. 


Polynomial Approximations 


The polynomial approximation method is fundamentally very simple. A limited part 
of a function is approximated by a polynomial of some order sufficient to obtain the desired 
accuracy. The polynomial is generally a series of the form: 


n 


P(n, x) = Dp [alilxi, (1) 
i=0 


where x is the independent variable, n the polynomial order (a fixed integer), and a[i] 
is a set of n+1 fixed coefficients. 


The desired function, say f(x), is then approximated by a particular P(n, x) such that: 
f(x) = P(n, x) + e(x), xl < x < xu, (2) 


where x1 and xu are the limits of the domain of x, and e(x) or e(x)/f(x) is the error func- 
tion which has been usually minimized in the min-max (equi-ripple) sense. This is done 
by selecting an appropriate means of calculating the coefficients ali]. 


Various techniques and schemes are used in.the selection of: 


e the approximation interval, 

e transformations on the function, 
e selection of the polynomial form, 
e error minimization criteria, and 
e calculation of the coefficients. 


See Hastings [2] for an excellent tutorial on this numerical methodology. All of the 
polynomial approximations used in here were obtained from the National Bureau of Stan- 
dards reference edited by Abramowitz and Stegun [3]. 


Non-Linear Equation Approximation 


The second method of approximation, using the solution of non-linear equations, 
is easier to understand. This method requires that a solution for the equation g(x) = 0 
be found. One means for solving this equation is by Newton-Raphson iteration. This can 
be understood by considering the Taylor series expansion for g(x): 


g(x + h) = g(x) + hg’(x) + r(x, h), - (3) 


where r(x, h) is the remainder of the series (which can be assumed to be small), and g’(x) 
is the derivative of the function g(x). Leaving off the remainder in (3) we get, in terms 
of incremental values of x, the approximation: 


g(xfi+1]) = g(xfi]) + {xfi+1]—x[i}g’(xfi)). (4) 
Solving for x[i+1] in (4) with g(x[i+1]) = 0 yields the approximation: 
x{it+1] = x[i] — g(xfi])/g’(xfi). (5) 


Thus, x[i+1] will converge to a solution of g(x) = 0. Convergence can be shown 
to be quadratic, i.e. the error in the approximation at each iteration is proportional to the 
square of the error in the previous iteration. Minimally, this requires a sufficiently close 
starting value for x[0] and the condition that |g’(x)| > 0 for all iterated values of x. 


Math Functions Details 


The approximation techniques can be applied to each of the classes of functions. 
The following sections describe the approximations as they are applied to each function. 


Inverse and Square Root Functions 


For the problem of computing good approximations to sqrt(c) (SQRT and SQRTX 
routines) and 1/c (FPINV and FPINVX routines), both g(x) and g’(x) must be derived 
and then use the iteration of equation (5). This is complicated by the restriction that divi- 
sion should be avoided since the TMS320C30 has no divide instructions. For the iteration 
to find the inverse of c, you can write: 


g(x(i]) = 1/xfi] — c = 0, (6) 


which is solved when 1/x = c or x = 1/c. Taking the derivative of (6) and substituting 
into (5) and simplifying gives us: 


x{i+1] = x[i](2 — cx[i}}, (7) 
which needs no division. 


Thus, (7) will converge to 1/c with the accuracy (in digits) for each iteration equal 
to twice that of the preceding one. Thus, if x[0] approximates 1/c to 3 bits of precision, 
only three iterations of (7) will yield about 24 = 3(23) bits of accuracy. 


A similar iteration from f(x) = x2 for sqrt(c) can be derived from the formulation: 

g(xli]) = xfi]? — ¢ = 0, (8) 
which is solved when x2 = c or x = sqrt(c). The solution for (8) leads to the classic square 
root formula: 

x{i+1] = 0.5{c/x[i] + x[i]}, (9) 


but this equation uses division. However, the iteration from f(x) = 1/x2 for 1/sqrt(c) can 
be shown to be: 


x{it+1] = x[i]{1.5 — c’x[i]2}, (10) 


where c’ = c/2 = 0.5c. Though (10) needs no division, the final desired result must be 
transformed by an extra multiplication by the input c because: 


sqrt(c) = c{1/sqrt(c)}. (11) 


Formula (10) will also converge, in the precision doubling fashion of the Newton- 
Raphson iteration, given a suitable close starting value for x[0] and the use of sufficiently 
accurate arithmetic. Note that the extended-precision version routines FPINVX and SQRTX 
both use an extra iteration (for a total of 4) to achieve the needed 32-bit accuracy for the 
40-bit format. 


The initial guess x[0], for the iterations of 1/sqrt(c) and 1/c, may be obtained using 
an interesting approximation. A TMS320C30 floating-point number c = (1 + m)2¢, where . 
0 < m < 1 and -127 < e < 127. The extra 1, added to the fractional mantissa m, 
is the implied bit. Then we can write the inverse of c as: 


I/e = 1/11 + m)2~-¢. (12) 
An excellent approximation for the inverse of the mantissa is: 
1/1 + m) = 1 —- m2, (13) 


which is exact at the end points: m = 0 and m = 1. Then the approximation for the 
reciprocal would be: 


l/c = (1 — m/2)2-¢. (14) 


It turns out that this approximation can be achieved in a single logical operation. 
If you compute the unlikely value of c’ = c XOR OFF7FFFFFFFh, you would comple- 
ment all bits in c except the sign bit. Including the implied bit and taking the effect of 
one’s complement arithmetic into account results in a final value of: 


ec’ = {1 + (1 — mj2-€ + D, (15) 
or the desired approximation: 
ce = (1 —m/)2-¢ = Ic. (16) 


c’ gives about 3 bits of precision, which is an excellent seed x[0] for the 1/c iteration. 
Using e/2, you have a start for the 1/sqrt(c) iteration as well. 


Sine and Cosine Functions 


The SIN, COS, SINX, and COSX (sine and cosine) routines all use the same basic 
approximation (section 4.3.98, p. 76 in [3]). The series is for sin(x)/x but is obviously 
transformed by multiplying by x. The polynomial of even terms then is of the form: 


5 
sin(x) = x )) {al2i}x2i} + xe(x), (16) 
i=0 


where |x| < Pi/2 and |xe(x)| < 2(10-9). Instead of using another power series for cos(x), 
you can use the fact that: 
cos(x) = sin(x + Pi/2). (17) 


The series given by (16) is only accurate in the 1st and 4th quadrants, i.e. Ix| < 
Pi/2. Sin(x) in the other two quadrants is found from: 


sin(x) = sin(Pi — x). (18) 


The case for x < 0 is expediently handled by using |x| for all calculations except 
for the final multiply by x in (16). 


Exponential Functions 
The EXP and EXPX (exponential) routines use an approximation (see Section 4.2.45, 
p. 71, in (3]). The expansion is of the form 


7 
exp(x) = )) {ali}xi} + e(x), (19) 
i=0 


where 0 < x < In(2) and |e(x)| < 2(10-!9). The series for 2¥ is found by substituting 
y = x/In(2) since: 


exp(x) = exp(In(2)y) = 2y. (20) 


The new expansion then becomes: 


7 
2y = ) [blilyi} + e(x), (21) 
i=0 
where b[i] = afi](In(2)i). See the coefficients in the EXP routine. 


Values of exp(x) for x outside the convergent range are found by two means. First 
for x < 0, note the relationship: 


exp(—x) = 1/exp(x), (22) 


which does require an inverse (see the FPINV and FPINVX routines). For y > 1, let 
y = n+ fwheren = 1,2,...and0 < f < 1. By substituting y in (20), you get 


exp(x) = 2n+f = (2f)(2n). (23) 
Natural Log Functions 
The LN and LNX (natural or base e logarithm) routines use the approximation from 


(3] (section 4.1.44, p. 69). The expansion comes in the form: 


8 
In(l + x) = )) {ali}xi} + e(x), (24) 
i=1 
where 0 < x < 1 and |e(x)| < 3(10-8). The expansion for In(y) can be used if the 
transformation y = x — | is applied. 


Values of In(x) for x outside the convergent range are found in the following way. 
First, make the substitution x = f(2") for 1 < f < 2andn=0,1,... , and then write: 


log2(x) = log2(f2") = n + log2(f), (25) 
where log2(x) is the log base 2 of x. Using the relationship that log2(x) = In(x)/In(2), 
you get the equation 

In(x) = In(f) + nln(2). (26) 
Arctangent Functions 


The ATAN and ATANX (arc or inverse tangent) routines use the approximation 
from section 4.4.49, p. 81 in [3]. The series with only even terms for atan(x)/x is trans- 
formed to 


8 
atan(x) = x Da {a[2i}x2i} + xe(x), (27) 
i=0 


where —1 < x < 1 and |xe(x)| < 2(10~8). Values for atan(x) for x outside the con- 
vergent range are obtained by noting the following identity: 


atan(x) = atan((x — 1)/(x + 1)) + Pi/4. (28) 


Using the bilinear transformation y = (x — 1)/(x + 1) assures, at the expense of 
a divide operation, that y < 1 for x > 1. The case for x < 0 is expediently handled 
by using |x| for all calculations except for the final multiply by x in (27). 


Divide and Multiply Functions 


The last group of routines in category I and II are those for the additional arithmetic 
functions FDIV and FDIVX (floating-point divides), and FMULTX (extended-precision 
floating-point multiply). The divide operation for the TMS320C30, a = b/c is done by 
calculating the reciprocal or inverse of the divisor c. Then you compute 


a = D(1/c). (29) 


For a normal-precision divide, FDIV finds 1/c by a call to FPINV. A subsequent 
normal TMS320C30 floating-point multiply of the rounded inverse provides a suitable 
quotient. For an extended-precision divide, FDIVX finds 1/c by a call to FPINVX. The 
inverse is then extended-precision multiplied by the dividend using FMULTX. 


The extended-precision floating-point multiply simulated by FMULTX is the key 
to the implementation of virtually all of the extended-precision functions. The extended 
multiply is achieved using the normal floating-point multiply of the TMS320C30. For two 
extended-precision numbers xa and xb, you can represent each as the sum of two floating- 
point numbers: xa = a + ea(2—24) and xb = b + eb(2—24). The quantities ea and eb 
are the one-byte extensions of xa and xb respectively. 


Thus the complete product xc = (xa)(xb) can be expanded and written as 
xc = (a)(b) + [(a)(eb) + (b)(ea)]2—24 + (ea)(eb)2 — 48. (30) 


The last term in (30) is always less than the 32-bit precision in the mantissa of the 
final result. Therefore, you need only to compute the first two terms in the product xc. 
Also, note that all the indicated products in (30) may be computed using a normal-precision 
native TMS320C30 multiply as long as the terms are collected in extended-precision 
registers. The additions are also done using the native TMS320C30 add as it is implemented 
in extended-precision. 


Integer Arithmetic Program Details 


Integer routines differ from the floating-point versions because they produce only 
integer results. If the computation can produce fractional values, then the fraction must 
be truncated to leave only the integer result. 


Integer Result Log Base 2 


The routine ILOG2 is a useful utility for computing integer value m of the log base 
2 of the integer n. The result is computed by successive multiplies by 2 (implemented 
as shifts by 1). The resulting relationship is n < 2m, such that if log2(n) is not an exact 
integer, m is rounded up to the next largest integer. This is useful as it allows the deter- 
mination of m from any value n > 0 (e.g. not a power of two) which might require the 
padding of additional values (zeros) for a radix 2 FFT. This program is very fast because 
of a delayed branch loop and internally requires only 4(m+1) cycles (cached) to do the 
calculation. 


Extended Precision Integer Multiply 


The IMULT routine is a modified version of the program EXTMPY in the 
TMS320C3x User’s Guide [1]. It has been modified and slightly speeded up. The negation 
of the final 64-bit product is done in two instructions by direct two’s complement nega- 
tion rather than by using one’s complement to simulate the same result. The product is 
computed by breaking the multiplier and multiplicand up into two 16 bit integers each. 
Thus the full product ¢ of the numbers a = au(2!6) + al, and b = au(2!6) + bl is 


c = (au)(bu)232 + [(au)(bl) + (bu)(al)]2!6 + (all), (31) 


where the powers of two indicated are accomplished by shifts. Note that each product 
in (31) must be represented as a 32-bit integer. The adds in the sum must be done with 
care to facilitate the carry between the two final 32-bit components of the product. 


Integer Divide 


The IDIV routine is a modified version of the program DIVI in the TMS320C3x 
User’s Guide [1]. It has been modified to return the absolute value of the remainder of 
the integer division. The remainder was originally computed, but was discarded during 
the extraction process for the quotient. A few more instructions allow the extraction of 
both the quotient and remainder from the result of the SUBC process. The program IDIV 
may be used for the computation of the modulo function. The output of IDIV is the pair 
{q, |r|} = a/b, with the property: 


0 <r = (a modulo b) < a, (32) 


fora > O and b > 0. The complete relationship is, by definition, a = bq + r, for positive 
a and b. 


Vector Utility Routines 


Vector utilities are functions which operate on arrays of numbers. Some utilities, 
like dot products and convolutions, are simple. Other utilities, like those presented here, 
are more involved. 


Complex and Complex Conjugate Array Multiplies 


The array routine *CORMULT computes the point-by-point complex conjugate 
multiply of two complex arrays. If the arrays are cl and c2, and are of length n, then: 


cl[{k] < cl{kJconj(c2[k]), k = 1,...,n, (33) 


where «- means replaces. Each complex array is assumed to be stored as two separate 
arrays, i.e. {cl} = {x1, y1} and {c2} = {x2, y2}. In cartesian complex representation, (33) 
becomes 


(x1 + iyl) — (x1 + iyl)(x2 — iy?2), (34) 


where i represents the imaginary constant sqrt(—1). Separating the real and imaginary 
parts, we have: 


xl + xlx2 + yly2, yl <— ylx2 — y2xl (35) 


This operation can be used for the frequency domain correlation of two FFTs to imple- 
ment time domain correlation. 


On the other hand, the array routine *CONMULT computes the point-by-point com- 
plex multiply of two complex arrays. If the arrays are cl and c2, and are each of length 
n, then 


cl[k] < cl{k](c2[k]), k = 1,...,n, (36) 
In cartesian complex representation, (36) becomes 

(xl + iyl) <— (x1 + iyl)(x2 + iy2). (37) 
Separating the real and imaginary parts results in 

xl © xlx2 — yly2, yl <— ylx2 + y2x1. (38) 


This operation can be used for the frequency domain convolution of two FFTs to imple- 
ment digital filtering. 


Complex Array Bit Reversal 


The array routine *CBITREV executes an in-place bit reverse permutation on two 
arrays simultaneously. This operation is generally used for index scrambling before a DIT 
FFT (decimation in time, see CIFFT2), or after a DIF FFT (decimation in frequency, 
see CFFFT2) for index unscrambling. Therefore, *CBITREV is useful in permuting com- 
plex arrays stored as two separate arrays which are associated with radix 2 FFTs. The 
program uses the bit reverse indexing feature of the TMS320C30 to achieve this function. 
The loop in *CBITREV is nearly as efficient in permuting two arrays together as per- 
muting one array alone. This is due to the use of parallel load and store instructions and 
a delayed (single cycle) conditional branch. 


Floating Point Conversions 


The array routines *FMIEEE and *TOIEEE are vectorized versions of their original 
scalar counterparts FMIEEE and TOIEEE. Both routines do fast conversions from or 
to IEEE format by avoiding dealing with special rare cases. Also, both programs convert 
the numbers in the arrays in-place which destroys the original data. These array versions 
of the format conversion routines are much faster than calling the scalar version routines 
in a special loop. These routines also have their own internal, shared constant table for 
conversions. 


Vector Primitives 


The array routines *VECMULT, *CONMOV, and *VECMOV are a useful suite 
of efficient programs for simple array operations. The first routine, *VECMULT, per- 
forms the simple operation x[k] <— x[{k]c which is a scalar-vector multiply useful in uniform- 
ly scaling an array by a constant c. You can use this for scaling arrays after an inverse 
FFT by choosing c = 1/n. The next routine, *CONMOV, performs the operation 
x(k] « c which is useful in filling or initializing any portion of an array to a single cons- 
tant c. The last routine, *VECMOV performs the simple operation x[k] «+ y[k], an array 
move, and is, therefore, generally useful. 


FFT Routines 


This category contains the two complementary radix 2 complex FFT programs 
CFFFT2 and CIFFT2. These programs differ from previously available TMS320C30 FFT 
programs in that they operate on complex arrays which are stored as two separate and 
independent real arrays. Both routines do the FFTs in-place and do no index permutations 
or constant scaling (multiplication). Also these programs require only a 3/4 cycle exter- 
nal, pre-computed sine table. As with previous FFT programs, these, too, have a special 
multiply-less butterfly loop for the occurrence of unity twiddle or complex rotation factors. 


The routine CFFFT2 is a DIF radix 2 complex forward FFT program and thus 
assumes a normally indexed pair of input arrays. The output array is bit-reverse permuted 
and normally must be unscrambled to be of any use (see *CBITREV). The routine CIFFT2 
is a DIT radix 2 inverse FFT program and thus assumes a bit-reverse indexed pair of 
input arrays. A normally indexed complex frequency spectrum must be bit-reverse scrambl- 
ed before using CIFFT2 (again, see *CBITREV). On the other hand, the output from 
this inverse FFT is in normal indexed order, but lacks the traditional scaling by the factor 
of 1/n. Therefore, back-to-back calls of CFFFT2 and CIFFT2 will return the original 
complex array (in proper order) but multiplied by a factor of n. Consult the handbook 
by Burrus and Parks [4] for additional FFT algorithm details. 


Linear Algebra Routines 


The routines *SOLUTN and *SOLUTNX are the normal- and extended-precision 
implementations of the algorithm for solving simultaneous linear equations. This algorithm 
is the modified Gauss-Jordan elimination without (off diagonal) pivoting. This is a simple 
algorithm which is intended for use with well-conditioned systems of dense linear equa- 
tions of moderate size. Well conditioned means that the system of linear equations is linearly 
independent or non-singular. This subject and further algorithm details are to be found 
in chapter 2 of [5] by Press et al, or any other book on the numerical techniques of linear 
algebra. This algorithm is suitable for a wide range of problems requiring the solution 
of a system of linear equations, e.g. exact or least squares polynomial fitting. 


A simple system of linear equations has the form: 


A[1, 1)x{1] + Af[1, 2]x[2] +... + Af], nJx{n] = y[1], (39) 
A(2, 1)x[1] + A[2, 2]x[2] +... + A[2, n]x[n] = y[2], 


An, 1Ux{1] + A[n, 2Ix(2] +... + Aln, mixin = yin]. 


Symbolically, you may write A = Ali, j] as the n x n matrix of coefficients, and 
x = x[i] as the unknown independent variable (column) vector, and y = y[j] as the depen- 
dent variable (row) vector. Thus (39) can be written in short hand form as Ax = y or 
Ax — y = 0, where the multiplication indicated is a matrix-vector multiply. The fun- 
damental problem in linear algebra, then, is to find the solution vector x. In fact, you 
may desire to find the m different solutions to m sets of linear equations which share the 
same coefficient matrix A, i.e. Ax(k] = y[k], fork = 1,...,m. 


You can solve the general problem just stated by using *SOLUTN, or with more 
accuracy with *SOLUTNX. This is done by constructing a tableau B (table of coefficients) 
which is simply the coefficient matrix A (in row major storage format) with the negative 
of the y vector(s) appended (:) as m extra columns to A. Thus you would have B = A 
: —y, as your problem, where B is a n by n+m matrix and typically m = 1. Thus, for 
the common case of m = 1, the input array B can be written as: 


A[i, 1], Afi, 2], ..., Afl, nl, —y(1], (40) 
A[2, 1], A[2, 2], ..., A[2, n], —y([2], 
Ain, 1), Aln, 0] are re n], Lint 


After the *SOLUTN routine is executed, the matrix C = A’ : x appears, where 
the column(s) beyond the original coefficients A (the y[k] vectors) have been replaced 
by the solution vector(s) x[k]. The new matrix A’ is a partially computed version of the 
inverse of the matrix A. The complete inverse of A, which is normally computed by the 
standard Gauss-Jordan scheme, is rarely needed. Therefore, a faster modified algorithm 
has been used which does about half the work. 


This simple method used for solving systems of linear equations has two restrictions. 


1. As the pivoting operation (exchange of x and y variables) always starts with 
A[1, 1] and proceeds down the diagonal, A[1, 1] must be non-zero. This is 
because, in the exchange process, you must divide by the pivot element. A zero 
coefficient at A[1, 1] may be moved by reordering the variable indices by ap- 
propriately swapping rows and columns in A and in y. 


2. The maximum absolute value of the elements in A must be approximately uni- 
ty. This is necessary to assure that no pivot element is encountered which is 
smaller in magnitude than 10-8 for *SOLUTN, and 10-10 for *SOLUTNX. 
This restriction monitors the system condition and assures an adequately ac- 
curate solution, but the final solution should always be verified by substitu- 
tion. This is done by inspecting the elements of the error vector e = Ax — 
y computed by using the solution x, and the original A and y. 


Summary 


This report presented a set of routines that can be used in digital signal processing 


applications. The appendix contains the source code of these routines. This source code 
can also be obtained from the Texas Instruments Electronic Bulletin Board (713) 274-2323. 
If there are comments or corrections, please contact the author of this report: 


(1) 


(2) 


(3) 


(4) 


(5) 


Mr. Gary Sitton 
Gas Light Software 
5211 Yarwell 
Houston, TX 77096 
Tel (713) 729-1257 
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